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We evaluate imaginary time density-density correlation functions for a two-dimeirsioiral homo¬ 
geneous electron gas using the phaseless auxiliary field quantum Monte Carlo method. We show 
that such methodology, once equipped with suitable numerical stabilization techniques necessary 
to deal with exponentials, products and inversions of large matrices, gives access to the calcula¬ 
tion of imaginary time correlation functions for medium-sized systems; we present simulations of a 
number up to 42 correlated fermions in the continuum, using up to 300 plane waves as basis set 
elements. We discuss the numerical stabilization techniques and the computational complexity of 
the methodology. We perform the inverse Laplace transform of the obtained density-density corre¬ 
lation functions, assessing the ability of the phaseless auxiliary field quantum Monte Carlo method 
to evaluate dynamical properties of many-fermion systems. 


I. INTRODUCTION 

The homogeneous electron gas (HEG) is one of the 
most widely studied systems in condensed matter physics 
It represents a model of recognized importance, 
which offers the opportunity to explore the quantum be¬ 
havior of many-body systems on a fundamental basis and 
provides a ground test for several quantum chemistry [ 3 , 
many-body 3 and quantum Monte Carlo (QMC) 

[m methodologies. Furthermore, recent years have wit¬ 
nessed the realization of increasingly high-quality two- 
dimensional (2D) HEGs in devices of considerable exper¬ 
imental interest such as quantum-well structures [H, 01 
and field-effect transistors M- 

The accuracy of QMC calculations for the HEG is 
unavoidably limited by the well-known sign problem 
[il[3, arising from the antisymmetry of many-fermion 
wavefunctions. The vast majority of QMC simulations 
of many-fermion systems circumvent the sign pro blem 
relying on the Fixed-Node (FN) approximation [TtI. [l8l| . 
Methodologies based on the FN approximation provide 
very accurate estimations of ground state properties such 
as the kinetic and potential energy, and the static struc¬ 
ture factor. On the other hand, as the extension of the 
FN approximation to the manifold of excited states is less 
understood and established [I1113: the study of dynam¬ 
ical properties of many-fermion systems is a very active 
and challenging research field [l^ [2ll - l^ . 

In a recent work [T^, performing an extensive study 
of exactly solvable few-fermion Hamiltonians, we have 
shown that the phaseless Auxiliary Fields Quantum 
Monte Carlo (AFQMC) [l^, method can become 

an important tool for the calculation of imaginary time 
correlation functions. 

Motivated by this result, in the present work we apply 
the phaseless AFQMC method to the 2D HEG. In partic¬ 


ular, we focus on the high-density regime {ts < 2), since 
our previous study has revealed that the computational 
cost of the algorithm increases severely with rg, making 
a study at high Vg hardly practicable. This is due to 
the fact that, increasing Vg, the number of plane waves 
required to reach convergence in the basis set size be¬ 
comes larger, making cumbersome to perform the linear 
algebra operations required by the methodology. From 
a more physical point of view, the stronger correlations 
give rise to a more pronounced curvature in the wave 
function, which is accurately reproduced by a large num¬ 
ber of Fourier coefhcients. Nevertheless, the high-density 
regime is extremely interesting as the presence of the in¬ 
teraction leads to the emergence of important correlation 
effects, enhanced by the low dimensionality. We evaluate 
density-density correlation functions in imaginary time 
F[q,T) and we perform their inverse Laplace transform 
to extract information about the excitations of the sys¬ 
tem. We introduce and describe a method for stabilizing 
the calculation of imaginary time correlation functions in 
AFQMC, and present numerical tests demonstrating its 
accuracy. We finally assess the accuracy of the calcula¬ 
tions comparing AFQMC results with predictions within 
the random phase approximation (RPA) for hnite sys¬ 
tems im. We also compare AFQMC estimates with 
the results of Fixed-Node calculations 0113 , performed 
with a nodal structure encompassing optimized rational 
backflow correlations 0,0,01 • 


The paper is organized as follows: the phaseless 
AFQMC method is briefly reviewed in Section [III the 
results of the study are discussed in Section iHll and con¬ 
clusions are drawn in the last Section CYl 
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II. METHODOLOGY 


state I'I't) has non-zero overlap with |$o) the relation: 


A. The Model 


|$o) OC lim (5) 

r—voo 


The 2D HEG is a system of charged spin-i fermions in¬ 
teracting with the Coulomb potential and immersed in a 
uniform positively-charged background. For the purpose 
of studying the 2D HEG we simulate a system of N par¬ 
ticles moving inside a square region TZ of surface O = , 

employing periodic boundary conditions (PBC) at the 
boundaries of the simulation domain, in conjunction with 
an Ewald summation procedure . In the present work, 
energies are measured in Hartree units Ena, and lengths 
in Bohr radii as- The Hamiltonian of the system reads, 
in such units: 


A .'t - ^ .'t 


k+qa^p-q<;^P‘i^ka 


kcr 


q^O k(T 

pq 


where spin-definite plane waves: 

^ikr 


( 1 ) 




2^fc,Z^. = ±- (2) 


with r G 7^, w = 0,1 are used as a basis for the single¬ 
particle Hilbert space. The ground-state energy per par¬ 
ticle of the system is obtained adding, to the mean value 
of (HJ), the corrective constant term: 
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erfc(y^|n|) 

^ |n| 

-^^^2 I I 


n^O 


-3.900265 
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( 3 ) 


holds, eo being the ground state energy, that can be esti¬ 
mated adaptively following a common procedure in Diffu¬ 
sion Monte Carlo (DMC) calculations [l^. QMC meth¬ 
ods rely on the observation that the deterministic evolu¬ 
tion can be mapped onto suitable stochastic processes 
and solved by randomly sampling appropriate probabil¬ 
ity distributions. Determinantal QMC methods, such as 
the phaseless AFQMC, use a Slater determinant as trial 
state I'I't), typically the Hartree-Fock state, and map dH) 
onto a stochastic process in the abstract manifold 'S{N) 
of A^-particle Slater determinants. This association is 
accomplished by a discretization of the imaginary time 

propagator 

I^-t) n G N,5r = - 

( 6 ) 

and by a combined use of the Trotter-Suzuki decompo¬ 
sition MM , of the Hubbard-Stratonovich transforma¬ 
tion [3lLl45ll4^ and of an importance sampling technique 

MM on the propagator e , The result is: 




Giv-mr) 


where dg{ri) is a multidimensional standard normal prob¬ 
ability distribution, G(r]) is a product of exponentials of 
one-body operators and: 

2n[r7,4]=e-^-^'«(vI/T|G(r7-OI^T) (7) 


arising from the Ewald summation procedure employed 


4I| . The Hamiltonian ([T]) can be parametrized in terms 


of the dimensionless Seitz radius defined by: 


n 

N 


1 

n 


( 4 ) 


where n is the density of the system and qb the Bohr 
radius. This parametrization shows that the matrix ele¬ 
ments of the kinetic energy roughly scale as |fcp ~ rj^, 
and those of the potential energy as 1/H|q| ~ Thus, 
for increasing Seitz radius, the interaction part of H plays 
a more and more relevant role. 


B. The Phaseless AFQMC 


is a weight function depending on a complex-valued pa¬ 
rameter which is chosen to minimize fluctuations in 
2 n[r7,|] to first order in St. Equation ([7]) illustrates 
the mechanism responsible for the appearence of the sign 
problem in the framework of AFQMC: when the overlap 
between G{r] — o-nd the trial state vanishes mas¬ 

sive fluctuations occur in ©• In the method conceived 
by S. Zhang, the exact com plex -valued weight function 
appearing in © is replaced |l9l . [30l| by the approximate 
form: 


2 H[r7,|] ~ x max(0,cos(A6»)) 


where eiod'^) = Re 
tional, and: 


(I'tI'I') 


( 8 ) 

is the local energy func- 


To address the calculation of static and dynamical 
properties of the 2DHE^ we rely on the phaseless 
AFQMC method [l^, , that moves from observa¬ 

tion that the imaginary time propagator e~'^^ acts as a 
projector onto the ground state |4)o) of the system in the 
limit of large imaginary time. Therefore, as long as a trial 


A9 = Im 


log 


{^T\G{rj-0\^)) 

(4't|4') 


(9) 


The first factor corresponds to the real local energy ap¬ 
proximation, which turns © into a real quantity, avoid¬ 
ing phase problems rising from complex weights; the real 
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local energy approximation is implemented neglecting 
some fluctuations in the auxiliary fields [l^ . The second 
factor, together with the introduction of the shift param¬ 
eters, has been argued in 0,1111 to keep the overlap be¬ 
tween the determinants involved in the random walk and 
the trial determinant far from zero. In fact, the angle A0 
corresponds to the flip in the phase of a determinant dur¬ 
ing a step of the random walk: the term max(0, cos(A0)) 
is meant to suppress determinants whose phase under¬ 
goes an abrupt change, under the assumption Pill 
that such behaviour indicates the vanishing of the over¬ 
lap with the trial state. 

C. Hubbard-Stratonovich Transformation in the 
plane-wave Basis Set 

In general, the structure of G{ri) is specified through a 
procedure that might result lengthy and computationally 
expensive 0, [ll| . When spin-definite plane waves are 
used as a basis for the single-particle Hilbert space, a 
remarkable simplification derived in details in Appendix 
occurs in its calculation and leads to the following 
result: 



Figure 1. (color online) pictorial representation of the ran¬ 
dom walk in the manifold of A-particle Slater determinants 
D(A) (lavender surface). The figure points out that the imag¬ 
inary time propagator drives a Slater determinant 

I^'t) away from D{N), while the one-body propagators G{rj) 
preserve T){N). This permits to retrieve the analytically in¬ 
tractable state I’Ft) as a stochastic linear combi¬ 

nation of Slater determinants G(»7o™^)l^r) according to (1151) . 


D. Numeric Implementation 


( 10 ) 

with: 



— ^ (^o)fc 

k(7 




( 11 ) 


and, denoting pg the Fourier component of the local den¬ 
sity: 


Ai(q) = 


StT Pq P—q 


H|q| 


M{q) = 


27r ipq ^P—q 


H|q| 


The operators (fl^ will be henceforth written as: 


^{q) = “ 




( 12 ) 


(13) 


kpa 


with: 


iM<i))kp = 
i^2iq))kp = \l 


27r 5k.p-\-q ”b ^k.p- 


H|q| 


27r 1^k,p-\-q '^^k, 


(14) 


p-q 


n\q\ 


We remind that formulae (iTni),(ini) and (fT^ result from 
an exact calculation, immediately generalizable to all ra¬ 
dial two-body interaction potentials, and to all spatial 
dimensionalities. 


The observations outlined above give rise to a polyno- 
mially complex algorithm for numerically sampling ((5|) . 
a pictorial representation of which is given in Fig. [T] 
Several Slater determinants {|^o™^)}^=i’ called walk¬ 
ers, are initialized to the Hartree-Fock ground state, a 
filled Fermi sphere in the case of translationally invariant 
systems such as the 2D HEG, and given initial weights 
equal to 1. 

Subsequently, each walker is let evolve under the action 
of the operators (fTUll and its weight is updated through 
multiplication by dS]). An estimate for the ground state of 
the system is provided by the following stochastic linear 
combination of Slater determinants: 




|$o) ^ 




(w) 




w) 




W — 1 




(15) 


Since numeric calculations can be carried out on finitely- 
generated Hilbert spaces only, the numeric implemen¬ 
tation of the phaseless AFQMC algorithm requires the 
single-particle Hilbert space basis ([2]) of the system to be 
truncated, i.e. only the M lowest-energy plane-waves to 
be retained. 


E. Dynamical Correlation Functions 

The formalism outlined in III B[ Hi PI enables the calcu¬ 
lation of ground state properties, and also of the imagi¬ 
nary time correlation function (ITCF): 

= {^o\Ae-<^-^°)B\^o) (16) 
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of two single-particle many-body operators A, B. CH) is 
a purely mathematical function, related to the dynamical 
or energy-resolved structure factor: 

/ dt—($o|4l(t)i3|$o) (17) 

^ IR 

of A, B, a quantity appearing in linear response the¬ 
ory and providing precious information on the time- 
dependent response of the system to external fields. Dy¬ 
namical structure factors and ITCFs are related to each 
other, as revealed by their Lehmann representation, by a 
Laplace transform [8|. 


Within the AFQMC formalism, the issue of comput¬ 
ing ITCFs is complicated by the circumstance that the 
one-body operators A, B do not map Slater determi¬ 
nants onto Slater determinants, but on rather compli¬ 
cated states. Nevertheless, making use of the canoni¬ 
cal anticommutation relations between fermionic creation 
and destruction operators it is possible to show [l^ that: 

= J dgiri)Biri)G{rj) (18) 

where B{ri) is a suitable one-particle operator. In the 

case of the 2D HEG, it reads: 

(19) 

hpa 

where: 

B{ig) = (20) 

is defined through: 

i'B{'n))k = g-^(Wo)p 

^ \ ) kp 

( 21 ) 

By application of (fT^ and of the backpropagation tech¬ 
nique [ll, [H, [13, it is possible to express the ITCF 
Ba b ('^) mean value of a random variable over the 
random path followed by the walkers in the manifold of 
Slater determinants. 


Further details of this calculation procedure are re¬ 
ported in [T3 . For the purpose of the present work, it is 
sufficient to recall that the phaseless AFQMC estimator 


of ^(r) reads: 
F^^^irSr) ~ 


where: 




Z^w — 1 


w = l ijkl \^BP,m\^n / 

■^V'/n—1 ^n—1’ ■ • ■ 5 ’In—r 

^ \'ln-l Sn-li ■ ■ • I '/n-r ^n-r)lj 






( 22 ) 


Blj=G\Vn- in)... GHrjn+m-l - in+m-l)\^T) 

(23) 


and: 


V{ri 

=V{ri 


(w) 

n—1 


(w) 

n—1 


t(u') (“) _ cM \ _ 

^n—1’ • ■ • 5 Un—r ^n—r) 

1 f 1 

Sn—1/ ■ * ' \ In—r *^n—rJ 


(24) 


The estimator (1^^ is essentially a weighted average of 
suitably-constructed matrix elements; each walker w con¬ 
structs the matrix element and the weights 

involved in the weighted average (l22)l from two 
Slater determinants |4>^^) and two matrices 

> »7i-l - d-l) These objects are func¬ 
tions of the auxiliary fields configurations defining 
the random path followed by the walker in the manifold 
of Slater determinants, and their calculation is pictorially 
illustrated in Fig. O 

In the present work, we consider the imaginary time 
density-density correlation function: 




(25) 


which is the Laplace transform of the dynamical struc¬ 
ture factor S{q,uj). This quantity is notoriously related 
to the differential cross section of electromagnetic radi¬ 
ation scattering, and provides essential information for 
the quantitative description of excitations of the HEG, 
collective charge density flutuations, i.e. plasmons, end 
electron-hole excitations [1,13 • 


F. Numeric Stabilization 


In a previous work [T^] we have pointed out that, due to 
the presence of such matrix elements, the estimator (1^^ 
is negatively-conditioned by a form of numeric instability. 
The aim of the present section is to elucidate the origin 
of such phenomenon and to propose a method for stabi¬ 
lizing the calculation of ITCFs in AFQMC. The AFQMC 
estimator (1^ of the ITCF F^ ^(r) involves a weighted 
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Figure 2. (color online) pictorial representation of the phaseless AFQMC estimator for F^g{T), equation (I22II : g{T) is 

computed at r = rhr with r = 2, with n = 5 propagation steps and m = 3 backpropagation steps. The matrix 2? appearing 
in (|24l) is computed between the time steps n — r and n (at which is computed); the determinant |^'^^) is computed 

between the time steps n and n + m, and the weights and are computed at the time steps m + n — r and m + n. 


average, over the random paths followed by the walk¬ 
ers employed in the simulation, of a quantity in which the 

matrix elements of , vItX ~ ^i-r) 

of its inverse appear. In the reminder of the present sec¬ 
tion, these matrices will be referred to as V and 'D~^ for 
brevity. The matrices "D and 'D~^ need to be computed 
numerically, respectively as product of r matrices and 
inverse oiD. It is well-known that the numerical compu- 
tation of T> and 'D~^ introduces rounding-off errors 
which accumulate as r increases with detrimental impact 
on the results of the computation M- 

Rounding-off errors are particularly severe when the 
oo-norm condition number: 

KiV) = \\V\\ao\\V-^\\oo (26) 

of the matrix 2?, in which ||A||oo = max^- \ Aij\ denotes 
the oo-norm on the space of M x M complex-valued ma¬ 
trices, is large. For the systems under study, we observe 
a condition number roughly increasing as k(T>) ~ C[ for 
some constant Ci. The rapid increase of k( 2?) indicates 
that the numeric matrix inversion I (2?) used to estimate 
2?~^ might be ill-conditioned, an intuition that can be 
confirmed by studying the /igure of merit: 

\\E\\^ = \\1-VI{V)\\^ (27) 

For small r, is comparable with the machine pre¬ 

cision e = 10“^®; it then increases as for some con¬ 
stant C 2 and eventually saturates around 1. In appendix 
a qualitative explanation of the power-law increase 
of is provided. The gradual corruption of data re¬ 

vealed by the increase of ||2ii||oo reflects, as illustrated 
in figure ([3]), on the quality of the AFQMC estimates 
of ITCFs, which combine the matrix elements of V and 
I iV) as prescribed by (1331) . We propose to mitigate the 
numeric instability of the ITCF estimator by perform¬ 
ing a Tikhonov regularization [48j| of the numeric inverse 


X (27). Practically, the SVD of V is computed: 

V = C/diag(cri ... (28) 

and X (27) is obtained as: 

1(27) =Fdiag(CTi... (29) 

where di = is defined by a regularization param¬ 

eter A. Large singular values ai ^ X are mapped to 
CTi ~ while small singular values < A are kept be¬ 
low the threshold Particular care must be taken in 
choosing the regularization parameter A, since for small 
A the Tikhonov regularization is clearly ineffective, while 
for large A it provokes a severe alteration in X (27). On 
the other hand, an intermediate value of A prevents small 
errors in 27, associated to small singular values (Ti, to be 
dramatically amplified by the numeric inversion. 


The effect of the Tikhonov regularization has been 
probed considering the model systems of 2 particles in¬ 
troduced in , for which exact numeric solution of the 
Hamiltonian eigenvalue problem is feasible, and thus the 
ITCFs is exactly known. In figure ([3]) we show the ef¬ 
fect of the Tikhonov regularization (1331) on the ITCFs. 
The results show the existence of a broad interval of A, 
comprising the machine precision e = 10“^®, for which 
the Tikhonov regularization mitigates the numeric insta¬ 
bility affecting the AFQMC estimator of ITCFs without 
introducing any appreciable bias besides that introduced 
by the real local energy and phaseless approximations. 
The figure displays, in the upper and lower panel respec¬ 
tively, the statistical errors of the AFQMC estimations 
and the discrepancies with respect to the exact results for 
three different values of A. It is evident that, as the imag¬ 
inary time becomes large, the effect of the regularization 
is very important. 
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Figure 3. (color online) Effect of the Tykhonoff regularization 
(I29II on an ITCFs relative to a system of = 2 electrons 
with M = 21 basis set elements. Upper panel: statistical 
uncertainty affecting the AFQMC estimate of F{q, r) with 
A = (lavender solid lines), A = (green dashed 

lines) and A = 0 (orange dotted lines). Lower panel: bias 
affecting the AFQMC estimate of F{q,r). 


G. Computational Cost 


The AFQMC estimator of ITCFs should join numeric 
stability and low computational cost. The aim of the 
present section is to show that the computational cost of 
(l2^ is 0{M^), M being the number of orbitals consti¬ 
tuting the single-particle basis. The contribution to 
(|2^ brought by a single walker of index w reads: 


Fw — ^ ^ ~A.rs Fkl (30) 

ijklrs 


where the abbreviation: 


(•)» = 




{w) I , |v[/(“)\ 
BP,m\ 1 ^" / 




(31) 


has been inserted. The generalized Wick’s theorem [l9|, 
Id^ implies that: 




(1321) is most conveniently expressed, introducing the defi¬ 
nition Qij = {a\aj)w and recalling canonical anticommu¬ 
tation relations, as: 

) tj ; = QrsQij T Qrj {_^is Qis) (33) 
Combining (1^ and (1^ yields: 

Fw = 'Y2 •^rsBklQrsQij'Fik'DJ'j^ 
ijklrs 

ri BkiQrj'Fik'Di-^ 

ijklr 

— GrjGisF>ikF>l'^ 

ijklrs 


Despite its cumbersome appearance, dSl can be ef¬ 
ficiently evaluated computing the intermediate tensors 
VB, AG^ and at the cost of 0{M^) operations, 

and subsequently computing F^ as: 

+ y;(pB).,(p-'sq,„A, (35) 

Ur 

Ur 

at the cost of 0{M^) more operations. The calcula¬ 
tion of Fu, further simplifies for operators A whose ma¬ 
trix elements read Aij = Aj Si a{j) for some function 
a : {1...M} —>■ {1... M}. The density fluctuation oper¬ 
ator pq = J2ka falls witlun such Category. 

The complexity 0{M^) is the best allowed by the 
phaseless AFQMC methodology: in fact, the calculation 
of ITCFs requires at least 0{M^) operations to accu¬ 
mulate the matrix F, and the contractions in (I34() do 
not compromise this favorable scaling with the number 
of single-particle orbitals. 


III. RESULTS 


In the present work we have simulated paramagnetic 
systems of A = 18,26,42 electrons at rg = 0.1, 0.5,1; we 
show also results for A = 18 particles at = 2. Our 
calculations qualify the phaseless AFQMC as a practical 
and useful methodology for the accurate evaluation of 
F(<7,t), for systems of A = 0(10^) correlated fermions 
in the continuum. The complexity scales as (M be¬ 
ing the number of basis sets elements), and the abso¬ 
lute statistical error of F{q,T) can be kept at the level 
10“^ — 2.5 10“^ with moderate computational resources 
even at values of r ~ ^/Ep for Vg = 0.1,0.5,1 and 
T ~ 2.b/Ep for Tg = 2, Ep = l/r^ being the Fermi 
energy. 
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The number N of electrons constituting the system 
is inferior to that typically involved in QMC ground 
state calculations of bulk fermionic systems, but com¬ 
parable to that used in the context of excited-states cal¬ 
culations through imaginary time correlation functions 
evaluated via configurational QMC methods reported in 
literature!^ ■ 

The imaginary time steps used in our calculations were 
6t = 0.003,0.004,0.006,0.008 at = 0.1,0.5,1,2 
respectively. For each simulation, the number of plane- 
waves constituting the single-particle Hilbert space has 
been raised up to M = 300 according to the number of 
particles and to the strength of the interaction. For all 
calculations, it was verified that decreasing the time step 
and increasing the number of plane-waves had a negligi¬ 
ble effect on the ground state energy. To obtain correct 
estimates of ITCFs, it is necessary to perform a suffi¬ 
ciently large number m of backpropagation steps. How¬ 
ever, it is well-known that raising m can result in 
an increase in variance, which severely limits the possi¬ 
bility of extracting physical information from the long 
imaginary-time tails of the ITCFs. We have used a num¬ 
ber of backpropagation steps in the range m = 200 — 600. 
When m = 600 has proved insufficient, to avoid the in¬ 
creases in variance mentioned above, AFQMC estimates 
have been extrapolated to the m —> oo limit (data ob¬ 
tained by extrapolation will be henceforth marked with 
an asterisk). 


A. Imaginary time correlation functions and 
excitation energies 

For all the values of N and r^, an AFQMC estimate 
of the ITCF (PSII is produced according to the procedure 
sketched in section Bl El The obtained F{q,T) is shown 
in the upper panel of Figures [H El [6] and [7l As it is 
well-known (5ll |. it is highly non-trivial to extract physi¬ 
cal information from ITCFs. In the case of the HEG, the 
finite size of the systems under study induces to expect 
contributions to F{q,T) coming from excited states of 
the system, which are obtained from the ground state by 
creation of particle-hole pairs. The dynamical structure 
factor S{q,uj), the inverse Laplace transform of F{q,T), 
is thus expected to display multiple peaks correspond¬ 
ing to the excitation energies. This picture is confirmed 
by RPA calculations for finite systems, reported in Ap¬ 
pendix o 

The presence of multiple peaks complicates the task 
of performing the analytic continuation providing an 
estimation of 5'(q,w). Therefore, since the number 
of peaks grows rapidly with the wave-vector modulus 
|q|, we have limited our attention to the wave-vectors 
qi = (27r/L)(l,0) and q2 = (27r/L) (1,1). Notice that 
\qi\/kF = 0.707,0.5,0.447 and \q 2 \/kF = 1,0.707,0.632 
for N = 18,26,42 respectively. Naturally, Uf = 
is the Fermi wave-vector. These low-momentum excita¬ 
tions are very interesting also from a physical point of 


view, in connection with the well-known collective plas- 
mon excitation of the HEG. Therefore, we do not attempt 
to predict the full S{q,u!), but limit ourselves to extract 
the excitations energies and weights by fitting the evalu¬ 
ated ITCE to a sum: 

F(r)=^s.e-™^ (36) 

i=l 

of exponential functions with positive frequencies Wi and 
weights Si relying on the well-established Levenberg- 
Marquardt curve-fitting method . The number 
Nw < 3 of such frequencies and weights is that leading 
to the best fit. 

In EiguresIH El [6]and[3we show results relative to the 
simulation of paramagnetic systems at rg = 0.1,0.5,1,2 
respectively. Each figure contains data relative to the 
particle numbers N = 18, 26,42 and wave-vectors qi,q 2 - 
In the upper panel we show the estimated F{q, r), while 
in the middle and lower panels we show, for qi and q 2 
respectively, the obtained frequencies and weights, to¬ 
gether with the RPA results. The AEQMC estimations 
of the quantities Si^iUi are displayed as points with both 
horizontal and vertical statistical errors: the horizontal 
ones provide the uncertainties on the frequencies uJi of the 
excitations, while the vertical ones gives the error bars 
on the weights Si. The coordinates of the points give, 
naturally, the mean frequencies and weights. The statis¬ 
tical uncertainties on the quantities Si^uji are those yield 
by the fit procedure. The frequencies predicted by the 
RPA are represented as impulses with height equal to the 
corresponding weights. We see that, at Vg = 0.1, there 
is a close agreement between AEQMC and RPA predic¬ 
tions of both frequencies and spectral weights. Since it 
is known that, for small rg, RPA predictions are very ac¬ 
curate, such agreement provides a robust check for the 
reliability of AEQMC methodology in providing informa¬ 
tion about the manifold of excited states of the system. 
It is well known that, in the same situation, calcula¬ 
tions of F{q,T) based on the Eixed-Node approximation 
would give inaccurate results even if the nodal structure 
of the ground state wave function is known with very high 
accuracy. As rg increases, discrepancies appear between 
the two approaches. The presence of such discrepancies 
is naturally expected: none of the methodologies used 
in the present work is free from approximations. The 
approximations underlying RPA and AEQMC, in par¬ 
ticular, are quite different in nature and are expected 
to agree only in the limit of high density (very low rg). 
Incidentally, we observe that, due to the finite size of 
the systems investigated, we cannot provide quantitative 
predictions about the plasmonic mode, which becomes a 
well defined collective excitation in the thermodynamic 
limit. What we expect is that, as the system size be¬ 
comes large enough, of the several peaks that are present 
for the finite systems, one will acquire a dominant spec¬ 
tral weight, eventually becoming a well defined collective 
excitation. In order to further assess the quality of our 




results, in Tables mini and m we detail the comparison 
with the RPA and configurational Monte Carlo results 
for the ground state energy per particle, the static struc¬ 
ture factor S{q) = F{q, 0) and the static density response 
function: 

In the AFQMC calculations, x{q) is obtained using the 
parameters yield by the fitting procedure of the density- 
density correlations F{q,T). As far as ground state en¬ 
ergies per particle are concerned, as shown in Table [H 
at Tg = 0.1 the three methods give compatible results. 
As Tg increases, AFQMC estimates are always closer to 
FN than RPA, lying between them. The configurational 
Quantum Monte Carlo evaluation of the ground state per 
particle has been obtained via Diffusion Quantum Monte 
Carlo (DMC) calculation, with a nodal structure encom¬ 
passing backflow correlations, optimized by means of the 
Linear Method [s^, . It is well-known that FN calcu¬ 

lations with optimized nodal structures yield highly ac¬ 
curate estimates of the ground state energy, as confirmed 
by comparison with Full Configuration Interaction QMC 
calculations this result, therefore, confirms the 

great accuracy of the phaseless approximation . 


The configurational QMC evaluation of the static 
structure factor S{q) has been obtained via FN DMC 
calculations with the nodal structure described above. 
DMC estimates have been extrapolated [i^. Moreover, 
we have compared AFQMC estimations of the static den¬ 
sity response function x(q) with RPA and Fixed-Node es¬ 
timations. In principle, the Fixed-Node evaluation of the 
static density response function x{q) is highly non triv¬ 
ial, involving the manifold of excited states. However, 
it is well-known that this difficulty can be circumvented 
m extracting x{q) from the ground state energy E{vq) 
of a system subject to an external periodic potential of 
amplitude Vq in the Uq —>■ 0 limit. We observe that, in¬ 
creasing Tg above 0.1, the AFQMC predictions remain, 
in general, closer to the configurational Monte Carlo ones 
than to the RPA ones: this is a strong indication about 
the quality of AFQMC results, since the Monte Carlo 
calculations include correlations beyond the RPA level. 
This result is remarkable, since the AFQMC evaluation of 
x{q) is considerably influenced by the low-energy excita¬ 
tions which, if predicted unaccurately, can significantly 
bias the result. As rg further increases, however, the 
agreement decreases. We have verified that the number 
of plane-waves and the number of backpropagation steps 
are sufficiently large to extrapolate the results and to fil¬ 
ter the excited states contributions from the trial wave 
function. Hence, the origin of the discrepancies between 
the estimations yield by the three methodologies used in 
the present work has to be sought in the approximation 
schemes underlying them. 


V rg ^ (RPA) 

a (AF) 

a (FN) 

18 0.1 40.14 

26 0.1 45.84 

42 0.1 42.18 

40.14(2) 

45.82(1) 

42.18(1) 

40.13(1) 

45.81(1) 

42.17(1) 

18 0.5 0.5065 

26 0.5 0.7520 

42 0.5 0.6031 

0.5007(2) 0.5012(2) 
0.7360(2) 0.7326(8) 
0.6002(1) 0.5922(9) 

18 1.0 -0.2489 
26 1.0 -0.1847 
42 1.0 -0.2215 

-0.2562(1) -0.2580(1) 
-0.1921(1) -0.1961(1) 
-0.2283(1) -0.2309(2) 

18 2.0 -0.2661 

-0.2695(1) -0.2717(1) 


Table I. RPA (column 3), AFQMC (column 4) and FN-DMC 
(column 5) estimates of the ground state energy for various 
systems (parameters are listed in columns 1-3); energies are 
measured in Ena- The RPA ground state energy is calculated 
on the Gaskell trial wavefunction [s^ . 


N Ts |q| 

S{q) (RPA) S{q) (AF) S{q) (FN) 

18 0.1 8.355427 

0.3105 

0.314(2) 

0.319(4) 

18 0.1 11.81636 

0.5150 

0.525(4) 

0.521(4) 

26 0.1 6.952136 

0.3326 

0.342(2) 

0.343(4) 

26 0.1 9.831805 

0.3623 

0.367(6) 

0.370(5) 

42 0.1 5.469911 

0.2101 

0.212(7) 

0.217(4) 

42 0.1 7.735622 

0.3045 

0.310(6) 

0.306(5) 

18 0.5 1.671085 

0.2511 

0.258(1) 

0.266(4) 

18 0.5 2.363271 

0.4137 

0.440(3) 

0.448(5) 

26 0.5 1.390427 

0.2225 

0.254(3)* 

0.238(4) 

26 0.5 1.966361 

0.3009 

0.313(2) 

0.322(4) 

42 0.5 1.093982 

0.1533 

0.161(2) 

0.146(5) 

42 0.5 1.547124 

0.2366 

0.247(2) 

0.264(4) 

18 1.0 0.835543 

0.2098 

0.231(2) 

0.218(5) 

18 1.0 1.181636 

0.3451 

0.395(3) 

0.386(4) 

26 1.0 0.695214 

0.1746 

0.227(2)* 

0.192(5) 

26 1.0 0.983181 

0.2558 

0.289(2) 

0.281(4) 

42 1.0 0.546991 

0.1219 

0.141(1) 

0.126(5) 

42 1.0 0.773562 

0.1938 

0.219(2) 

0.208(5) 

18 2.0 0.417771 

0.1657 

0.172(2)* 

0.176(4) 

18 2.0 0.590818 

0.2732 

0.304(3)* 

0.305(4) 


Table II. RPA (column 4), AFQMC (column 5) and FN-DMC 
(column 6) estimates of the static structure factor S{q) for 
various systems and wave-vectors (parameters are listed in 
columns 1-3); wave-vectors are measured in AFQMC 

estimates marked with an asterisk are extrapolated. 


IV. CONCLUSIONS 

We have shown the possibility to provide accurate first 
principles calculations of imaginary time correlations for 
medium-sized fermionic systems in the continuum, us¬ 
ing the phaseless Auxiliary Fields Quantum Monte Carlo 
method. 
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Figure 4. (color online) Upper panel: imaginary time correlation functions of the density fluctuation operator pq for param¬ 
agnetic systems of = 18, 26 and 42 particles (left to right) at = 0.1, with transferred momenta q\ (green dashed lines) 
and q 2 (lavender solid lines). When not visible, errors are below the symbol size. Lines are only a guide for eyes. Central 
panel: dynamical structure factor for N = 18, 26 and 42 particles (left to right) with transferred momentum qi (RPA: orange 
impulses, AFQMC: green symbols). Lower panel: dynamical structure factor for N = 18,26 and 42 particles (left to right) 
with transferred momentum q 2 (RPA: orange impulses, AFQMC: lavender symbols). 


We have simulated a 2D homogeneous electron gas of 
up to = 42 electrons using a plane-waves basis set of up 
to M = 300 elements. We have shown that the density- 
density correlation function in imaginary time can be cal¬ 
culated via a polynomially complex algorithm with the 
favorable scaling 0{M^). In order to achieve a good ac¬ 
curacy level in the calculations, we propose stabilization 
procedures to deal with matrix inversion, which can be 
used in combination with well-established stabilization 
procedures for matrix exponentiation and multiplication 
in particular, we suggest a Tikhonov regulariza¬ 
tion that allows to maintain a good accuracy level even 
for imaginary time values of the order of 3/Ep. We have 
yielded also comparisons with predictions of the static 
structure factor and the static density response obtained 
via the RPA approximation and via Fixed-Node Quan¬ 


tum Monte Carlo calculations. 


At small Vs the AFQMC correctly reproduces the RPA 
results. At larger Vg on the other hand, it provides quan¬ 
titative estimates of the deviations from the RPA, as 
the comparison with FN calculations reveals. We be¬ 
lieve this is a relevant result for QMC simulations: it is 
known, in fact, that the widely employed Fixed-Node ap¬ 
proximation fails to properly sample the imaginary-time 
propagator, due to the imposition of the ground-state 
nodal structure to excited states. AFQMC, on the other 
hand, appears to provide a useful tool to explore, from 
first principles, the manifold of the excited states of a 
fermionic system. 



























10 


N Ts \q\ 

x{q) (RPA) x{q) (AF) xiq) (fn) 

18 0.1 8.355427 

0.00276 

0.0028(4) 

0.00287(1) 

18 0.1 11.81636 

0.00449 

0.0046(1) 

0.00469(1) 

26 0.1 6.952136 

0.00598 

0.0065(2) 

0.00653(4) 

26 0.1 9.831805 

0.00282 

0.0028(1) 

0.00287(4) 

42 0.1 5.469911 

0.00311 

0.0032(4) 

0.00325(4) 

42 0.1 7.735622 

0.00330 

0.0034(2) 

0.00335(2) 

18 0.5 1.671085 

0.04516 

0.048(1) 

0.0484(4) 

18 0.5 2.363271 

0.06992 

0.081(2) 

0.0788(4) 

26 0.5 1.390427 

0.06298 

0.085(6)* 

0.069(2) 

26 0.5 1.966361 

0.04827 

0.051(1) 

0.0504(4) 

42 0.5 1.093982 

0.04074 

0.043(5) 

0.042(2) 

42 0.5 1.547124 

0.04903 

0.051(3) 

0.049(2) 

18 1.0 0.835543 

0.12612 

0.152(3) 

0.143(2) 

18 1.0 1.181636 

0.18979 

0.301(6) 

0.226(1) 

26 1.0 0.695214 

0.14601 

0 . 22 ( 2 )* 

0.162(2) 

26 1.0 0.983181 

0.13872 

0.188(6) 

0.161(2) 

42 1.0 0.546991 

0.10212 

0.158(7) 

0.14(1) 

42 1.0 0.773562 

0.13014 

0.176(9) 

0.16(1) 

18 2.0 0.417771 

0.31451 

0.34(1)* 

0.374(4) 

18 2.0 0.590818 

0.46238 

0.89(1)* 

0.590(2) 


Table III. RPA (column 4), AFQMC (column 5) and FN- 
DMC (column 6 ) estimates of the compressibility x{q) for 
various systems and wave-vectors (parameters are listed in 
columns 1-3); wave-vectors are measured in and x(q) in 
AFQMC estimates marked with an asterisk are extrap¬ 
olated. 


and pq = J2ka is the density fluctuation opera¬ 

tor. Recalling the parity of ^ and the anticommutation 
relation: 


[Pq,P- 


(Pq+P-q)'^ {ipq-ip-qf 


-qj+ - 2 

one eventually finds: 

1 


H = 


with: 

and: 

A{q) = 


k(7 




(A3) 


i?o-t-- ^ (Ai(q)^-I-^ 2 ( 9 )^^ (A4) 


(A5) 


27r Pq + p^q 


Q|q| 


^ 2 ( 9 ) = 


27r ipq — ip. 


Q|q| 


(A 6 ) 

which, since p^q = pi, are hermitian operators. Ap¬ 
plying the Hubbard-Stratonovich transformation to the 
propagator of the Hamiltonian (111(1 yields: 


G[ri) = mqTi(q)+r;2,d2(q)g-^ilo 
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The distance: 

\\X{Vr)-V-^\\^ (Bl) 

between the actual inverse 2 ?“^ of TPr and its numeric 
estimate I {Vr) (IBIII is bounded by: 

\\X{Vr) - < M\\X{Vr) 


Appendix A: Hubbard-Stratonovich Transformation 
for the 2DHEG 

By a straightforward application of the canonical an- 
ticommutation relations, the Hamiltonian o can be ex¬ 
actly rewritten as: 

V—\ / |fcp ^ - * 

fcdT ^ ^ q 

(Al) 

where: 


with: 

\\Er\\oo = \\l-VrX{Vr)\\oo (B3) 

Equation (IB2I) holds for ||i?r||oo < ]g, and is therefore 
adequate to the description of ||i?r||oo for small r. It can 
be combined with the following estimate [s^ : 

]\/T3 

\\X[Vr)-V-^U^e\\V~^\\1— (B4) 

to yield: 


p.{k) 


— V 

2Q ^ 


p^fc 


27r 

\p-k\ 


(A2) 


.M^ 


M||E,||oo=:^ 


\n 


- 1||2 
00 


U{T^r)\l 




(B5) 
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In the case of AFQMC calculations, where Vr and ^ 
come from the product of r matrices: 

\\V-^U=Cl \\I{Vr)U=Cl (B6) 


where Ca and (74 are suitable constants, close to each 
other. Merging (IB4I) and (IB6I) leads to: 


ll^^rlloo 



(B7) 


which reduces to: 


llArIloo 




(B8) 


in the limit of small r. Since C/^<C\, C-^ and C/^ being 
close to each other, the estimate (IB8I) leads to a power- 
law increase of HAr-Hoo- 


Appendix C: RPA for Finite Homogeneous Systems 

The aim of this appendix is to provide a brief de¬ 
scription of the Random Phase Approximation (RPA) 
1,1113 for finite interacting systems and of the proce¬ 
dure leading to the excitation energies and weights, with 
which AFQMC results have been compared. 

The RPA can be regarded to [ 1 ] as a refinement of 
the well-known Tamm-Dancoff approximation IS 113 
(TDA), which has long been supporting the study of 
excitations in nuclear systems. The TDA relies on the 
assumptions that the ground state of the system is the 
Hartree-Fock determinant, and that excited states can be 
represented as superpositions of determinants obtained 
promoting a single particle above the Fermi surface. 
Within RPA, on the other hand, a better approxima¬ 
tion |$o) for the actual ground state of the interacting 
system is employed to build up an Ansatz for plasmonic 
wave functions. To this purpose, the distinction between 
spin-orbitals below and above the Fermi level is made 
explicit by writing: 


= 


if > kp 
if \k\ < kp 


(Cl) 


^kcr 


and the Hamiltonian 0 is consequently expressed as: 
H = f + V = 


k(7 k(7 

■f w X! ^qPqP-q 


- ) + 


(C2) 


q/O 


where tk = 4'q = fh® first sum goes over all 


wave-vectors k such that |fc| > kp, the second sum goes 
over all wave-vectors k such that |fc| < kp and the density 
fluctuation operator pq is approximated [13 by: 

Pq — ^ ^k+qa^ka + k^+qa^ka (^2) 

k(7 


where the first sum, describing forward scattering pro¬ 
cesses in which a particle is promoted above the Fermi 
level, goes over all wave-vectors k such that |fe| > kp and 
|fc -|- q| < kp, and the second sum, describing backward 
scattering processes in which a particle is brought back 
below the Fermi level, goes over all wave-vectors k such 
that |fc| < kp and |fc -|- q| > kp. The RPA Ansatz for 
plasmonic wavefunctions is: 


k(7 k(7 


is justified by the observation that the pair destruc¬ 
tion operator bk+qa-Cka- annihilates the Hartree-Fock de¬ 
terminant but not the actual ground state of the interact¬ 
ing system. The eigenvalues e such that H\^q) = el^q) 
are obtained recalling that the commutators between the 
Coulomb interaction and the pair creation and destruc¬ 
tion operators can be approximated [13 as: 


^k+qa^kcr’^] — ^ Pq 


and: 


{kk+qcr^ktT^^I — n Pq 


(C5) 


(C6) 


respectively. Now, since |<i)o) and |4>q) and eigenstates of 
H with eigenvalues cq and e = eo + Ae respectively, the 
following identity holds: 

0=(4>q|(e-H)4+^Xl‘fio) = 

= ^^{%\ci + qAj^o) - {%\{H,Aq;bij<^o) 


l^o) = 


from which: 

(‘fi<?|Cfc+qcr'^fc<T 
follows. Similarly: 

(41q I ^fc-j-qo-Cfecr 14*o) ~ 


Ae tk tk-\-q 




(C7) 


(C8) 


(C9) 


Equation (IC8I) and 
the secular equation: 


Ae -j- tk tk-{-q 

can be summed over fc, cr to yield 





12 


(^glPql^o) 


‘^4'q 


(^gl/3<?l^o) 


E 


\k\<kF 

\|fc+q|>fcF 


which, simplifying the matrix element ($g|/5g|$o) in both 
members, and applying the change of variables r = —k — 
q in the second sum, takes the form: 



2 


E 

I A:| 

|fc+q|>fcF 


1 

tfc tk-\-q ~t“ ^6 


1 

ik ik-\-q 


(Cll) 

where the term between square brakets is immediately 
identified with the real part of the 2D Lindhard func¬ 
tion xo(q, Ae) i- The coefficients Xk, Yk are determined 
substituting (IC4I1 in (IC8I) and (IC9p . and read: 


Xk 


Yk 


N 

tk tfc+q -j- Ae 
M 

tk tk+q "h Ae 


(C12) 


where A/” is a normalization constant. Notice that Xk is 
defined for |fc| < kp, |fc + q| > kp while Yk for |fc| > kp 
and \k + q\ < kp. The right-hand side of (IClip is a 
function /(Ae) = ^ xo( 9 i Ae), illustrated in Fig.[8l with 
the following properties: 


lim /(Ae) < 0 
lim /(Ae) = 0"^ 

Ae—>-+oo 


(CIS) 


and diverging in corrispondence to the particle-hole en¬ 
ergies — tk- As a consequence, there exists a root 


1 

Ae tk tk-\-q 


E 

\k\>kF 




Ae -t- tk tk-\-q 


{CIO) 


of the secular equation (ICllI) between all the poles of 
/(Ae) and another root above them. The excited state 
corresponding to this root has coefficients A^, Yk sharing 
the same sign, and is therefore a coherent superposition 
of particle-hole excitations describing a collective high- 
energy oscillation being precursive of the plasmon. The 
excited states orresponding to other roots of (ICllll have 
coefficients Xk,Yk with non-constant sign, and therefore 
take into account the persistence of non-interacting prop¬ 
erties in the spectrum of the electron gas, even in pres¬ 
ence of Coulomb interaction. 


We have seen that the RPA approximation yields an 
Ansatz for the energies eq^„ and wavefunctions |d)q_„) of 
excited states with definite momentum q, which results 
in the following approximation for the image of the RPA 
ground state through the density fluctuation operator pq\ 

Pql^o) =^|$q,n)($q,n|Pq|$o) (C14) 

n 

with: 

(‘bq,n|Pq|^o) = E Xk,n + E Yk,n (CIS) 

|fe|<A:ir |fe|)>A:jr 

\k-\-q\>kF \k-\-q\<kF 

and for the dynamical structure factor: 

*^(9’ ^ E (C16) 
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Figure 5. (color online) Upper panel: imaginary time correlation functions of the density fluctuation operator pq for param¬ 
agnetic systems of = 18, 26 and 42 particles (left to right) at = 0.5, with transferred momenta q\ (green dashed lines) 
and q 2 (lavender solid lines). When not visible, errors are below the symbol size. Lines are only a guide for eyes. Central 
panel: dynamical structure factor for N = 18, 26 and 42 particles (left to right) with transferred momentum qi (RPA: orange 
impulses, AFQMC: green symbols). Lower panel: dynamical structure factor for N = 18,26 and 42 particles (left to right) 
with transferred momentum q 2 (RPA: orange impulses, AFQMC: lavender symbols). 
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Figure 6. (color online) Upper panel: imaginary time correlation functions of the density fluctuation operator pq for param¬ 
agnetic systems of = 18,26 and 42 particles (left to right) at Ts = 1, with transferred momenta qi (green dashed lines) 
and q 2 (lavender solid lines). When not visible, errors are below the symbol size. Lines are only a guide for eyes. Central 
panel: dynamical structure factor for N = 18, 26 and 42 particles (left to right) with transferred momentum qi (RPA: orange 
impulses, AFQMC: green symbols). Lower panel: dynamical structure factor for N = 18,26 and 42 particles (left to right) 
with transferred momentum q 2 (RPA: orange impulses, AFQMC: lavender symbols). 
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Figure 7. (color online) Upper panel: imaginary time correla¬ 
tion functions of the density fluctuation operator pq for para¬ 
magnetic systems of A'’ = 18, 26 and 42 particles (left to right) 
at Vs = 2, with transferred momenta qi (green dashed lines) 
and <72 (lavender solid lines). When not visible, errors are 
below the symbol size. Lines are only a guide for eyes. Cen¬ 
tral panel: dynamical structure factor for N — 18, 26 and 42 
particles (left to right) with transferred momentum qi (RPA: 
orange impulses, AFQMC: green symbols). Lower panel: dy¬ 
namical structure factor for N = 18, 26 and 42 particles (left 
to right) with transferred momentum <72 (RPA: orange im¬ 
pulses, AFQMC: lavender symbols). 







17 



0.0 1.0 2.0 3.0 

Ae [Ha] 


Figure 8. (color online) RPA secular equation for a paramag¬ 
netic system of = 18 particles at = 1, and q 2 . The blue 
solid line is the right member /(Ae) of (ICllI) . and the orange 
dashed line is the constant function g(Ae) = 1; intersections 
between the two graphs are marked with red dots, the RPA 
eigenvalues are the abscissae of such intersections. 











